7-3 All Cause Mortality Day 28

Analyse of the all-cause mortality secondary outcome in ASCOT.

Authors

Affiliations

Rob Mahar

University of Melbourne

James Totterdell

University of Sydney

Preamble

Load packages

library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(lubridate)
library(patchwork)
library(ggdist)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))
bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))
color_scheme_set("red")
options(digits = 4)

Prepare dataset

devtools::load_all()
all_dat <- read_all_no_daily()

# Anticoagulation set (intention-to-treat, including missing)
acs_itt_dat <- all_dat %>% 
  filter_acs_itt() %>%
  transmute_model_cols_grp_aus_nz()

# Anticoagulation set (intention-to-treat, excluding missing)
acs_itt_nona_dat <- acs_itt_dat %>% 
  filter(!is.na(D28_death))  

# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
  filter_concurrent_intermediate()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>% 
  filter(!is.na(D28_death))

# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
  filter_concurrent_std_aspirin()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>% 
  filter(!is.na(D28_death))

# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
  filter_concurrent_therapeutic()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>% 
  filter(!is.na(D28_death))

Load the Stan models

model_full <- cmdstan_model(file.path("stan", "binary", "logistic_site_epoch.stan"))
model_site_only <- cmdstan_model(file.path("stan", "binary", "logistic_site.stan"))
model_epoch_only <- cmdstan_model(file.path("stan", "binary", "logistic_epoch.stan")) 
model_fixed_only <- cmdstan_model(file.path("stan", "binary", "logistic.stan"))

Helper functions for this analysis script

make_stan_data <- function(
    dat, 
    vars = NULL, 
    outcome = NULL, 
    beta_sd = NULL, 
    ctr = contr.orthonorm) {

  X <- make_X_design(dat, vars, ctr)
  nXassign <- sum(grepl("rand", colnames(X))) - 1
  beta_sd <- c(2.5, rep(1, nXassign), beta_sd)

  y <- dat[[outcome]]
  N <- dim(X)[1]
  K <- dim(X)[2]  
  
    epoch  <- dat$epoch
  M_epoch  <- max(dat$epoch)
    region <- dat[["ctry_num"]]
  M_region <- max(region)
    site <- dat[["site_num"]]
  M_site <- max(site)
  region_by_site <- region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)
  
  list(N = N, K = K, X = X, y = y,
       M_region = M_region, region = region,
       M_site = M_site, site = site,
       M_epoch = M_epoch, epoch = epoch,
       region_by_site = region_by_site,
       beta_sd = beta_sd)
}

make_stan_data_concurrent <- function(
    dat, 
    vars = NULL, 
    outcome = NULL, 
    beta_sd = NULL, 
    ctr = contr.orthonorm) {
  
  X <- make_X_design(dat, vars, ctr, includeA = FALSE)
  nXassign <- sum(grepl("rand", colnames(X))) - 1
  beta_sd <- c(2.5, rep(1, nXassign), beta_sd)
  y <- dat[[outcome]]
  N <- dim(X)[1]
  K <- dim(X)[2]  
  list(N = N, K = K, X = X, y = y, beta_sd = beta_sd)
}

make_mortality_table <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3", "C4"),
      labels = str_replace(intervention_labels()$CAssignment[-1], "<br>", " "))) %>%
    group_by(CAssignment) %>%
    summarise(
      Patients = n(),
      Known = sprintf(
        "%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
      Missing = sprintf(
        "%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
      `Died by day 28` = sprintf(
        "%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Patients = n(),
          Known = sprintf(
            "%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
          Missing = sprintf(
            "%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
          `Died by day 28` = sprintf(
            "%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(key, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(font_size = 10, latex_options = "HOLD_position")  %>%
    row_spec(nrow(tab), bold = TRUE)
}

Descriptive

Anticoagulation

Summary table

save_tex_table(make_mortality_table(
      all_dat %>%
      filter_acs_itt(),
      "latex"), 
      "outcomes/secondary/anticoagulation-summary-mortality")
make_mortality_table(
    all_dat %>%
      filter_acs_itt()
)

n (\%)	Patients	Known	Missing	Died by day 28
Low dose	610	596 (97.7)	14 (2.3)	19 (3.2)
Intermediate dose	613	603 (98.4)	10 (1.6)	15 (2.5)
Low dose with aspirin	283	281 (99.3)	2 (0.7)	10 (3.6)
Therapeutic dose	50	50 (100.0)	0 (0.0)	6 (12.0)
Overall	1556	1530 (98.3)	26 (1.7)	50 (3.3)

Summary of mortality to day 28 by anticoagulation intervention.

Age

Relationship age to outcome

agedat <- acs_itt_dat %>%
  dplyr::count(D28_death, AgeAtEntry) %>% 
  spread(D28_death, n, fill = 0) %>% 
  mutate(
    n = `0` + `1` + `<NA>`,
    p = `1` / (`1` + `0`))
agemod <- glm(
  cbind(`1`, `0`) ~ AgeAtEntry, 
  data = agedat, 
  family = binomial())
agedat <- agedat %>%
  mutate(
    ypred = predict(agemod, newdata = agedat, type = "response")
  )
p1 <- ggplot(agedat, aes(AgeAtEntry, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry")
p2 <- ggplot(agedat, aes(AgeAtEntry, p)) +
    geom_point() +
    geom_vline(xintercept = 60, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\ndied by day 28", x = "Age at entry")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-age.pdf"), p, height = 2, width = 6)
p

Figure 1: Relationship (logistic regression linear in age) between age at entry and the primary outcome.

Country

Relationship country to outcome

tdat <- all_dat %>%
  filter_acs_itt() %>%
  dplyr::count(Country = fct_infreq(PT_CountryName), D28_death) %>%
  group_by(Country) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(tdat, aes(Country, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Country of enrolment")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-country.pdf"), p, height = 2, width = 6)
p

Figure 2: Mortality by country.

Site

Relationship site to outcome

tdat <- all_dat %>%
  filter_acs_itt() %>%
  dplyr::count(
    Country_lab = Country,
    Site_lab = fct_infreq(Location),
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")),
    Site = PT_LocationName,
    D28_death) %>%
  group_by(Country, Site) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  ungroup()

p1 <- ggplot(tdat, aes(Site_lab, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(tdat, aes(Site_lab, p_1)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Site of enrolment") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p <- p1 / p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-country-site.pdf"), p, height = 4, width = 6.25)
p

Figure 3: Day 28 mortality by site within country.

Calendar Time

Relationship calendar date to outcome

caldat <- all_dat %>% 
  filter_acs_itt() %>%
  dplyr::count(D28_death, yr = year(RandDate), mth = month(RandDate)) %>% 
  spread(D28_death, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`),
         n = `1` + `0` + `<NA>`)
p1 <- ggplot(caldat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12) +
  ylim(0, 0.15)
p2 <- ggplot(caldat, aes(mth, n)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p <- p2 | p1
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-calendar-time.pdf"), p, height = 2, width = 6)
p

Figure 4: Relationship between calendar date and the primary outcome.

Modelling

The analysis for all cause mortality is specified equivalently to that for the primary outcome. It includes intervention as randomised, age category, country, site nested within country, epoch, and intervention eligibility. The main analysis restricts the population to only those participants randomised to the anticoagulation domain.

Primary model

Data and sampling

seed <- 59579
mdat <- make_stan_data(
    dat = acs_itt_nona_dat, 
    vars    = c("inelgc3", "agegte60", "ctry"),
    outcome = "D28_death",
    beta_sd = c(10, 2.5, 1, 1))
snk <- capture.output(
  mfit <- model_full$sample(
    data = mdat,
    seed = seed,
    refresh = 0,
    iter_warmup = 1000,
    iter_sampling = 2500,
    chains = 8,
    parallel_chains = 8,
    adapt_delta = 0.99
))
mfit$save_object(file.path("outputs", "models", "secondary", "7-3.rds"))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(mdat$X)
names(mdrws$gamma_epoch) <- acs_itt_nona_dat %>% dplyr::count(epoch, epoch_lab) %>% pull(epoch_lab)
names(mdrws$gamma_site) <- acs_itt_nona_dat %>% dplyr::count(site_num, site) %>% pull(site)
site_facet <- factor(mdat$region_by_site, labels = c("India", "Australia\nNew Zealand", "Nepal"))

mdrws$Acon <- attr(mdat$X, "contrasts")$randA %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))]
mdrws$Atrt <- mdrws$Acon[-1] - mdrws$Acon[1]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c(
  "Intermediate vs low" = exp(mdrws$Ctrt[1]),
  "Low with aspirin vs low" = exp(mdrws$Ctrt[2]),
  "Therapeutic vs low" = exp(mdrws$Ctrt[3]),
  "Intermediate vs low with aspirin" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]),
  "Intermediate vs therapeutic" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[3]),
  "Low with aspirin vs therapeutic" = exp(mdrws$Ctrt[2] - mdrws$Ctrt[3])
)
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic")),
  setNames(exp(mdrws$beta[7:8]), c("Ineligible aspirin", "Age \u2265 60")),
  setNames(exp(mdrws$beta[9:10]), c("Australia/New Zealand", "Nepal"))
)

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-3-primary-model-acs-itt-summary-table")
odds_ratio_summary_table(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate	0.84	(0.42, 1.70)	0.89 (0.33)	0.69
Low with aspirin	0.82	(0.35, 1.84)	0.89 (0.39)	0.69
Therapeutic	3.17	(0.96, 10.05)	3.77 (2.43)	0.03
Ineligible aspirin	4.73	(1.23, 16.96)	5.82 (4.21)	0.01
Age ≥ 60	2.08	(1.12, 3.86)	2.19 (0.72)	0.01
Australia/New Zealand	0.62	(0.13, 2.75)	0.82 (0.73)	0.73
Nepal	2.85	(0.57, 11.08)	3.61 (2.87)	0.09

Epoch and site terms

p <- plot_epoch_site_terms(
  exp(mdrws$gamma_epoch), 
  exp(mdrws$gamma_site), 
  site_facet
)
pth <- file.path(
  "outputs", "figures", "outcomes", "secondary", "7-3-primary-model-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Code

p <- plot_or_densities(mdrws$compare)
p

Posterior densities for treatment comparisons.

Diagnostics

Code

mfit$summary(variables = c("beta", "gamma_site", "gamma_epoch"))

# A tibble: 48 × 10
   variable    mean  median    sd   mad     q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>  <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -4.73   -4.68   0.765 0.754 -6.05  -3.56    1.00   11020.   13326.
 2 beta[2]   0.975   0.977  0.483 0.482  0.189  1.78    1.00   24918.   15667.
 3 beta[3]  -0.222  -0.220  0.282 0.281 -0.689  0.239   1.00   26301.   16533.
 4 beta[4]  -0.499  -0.501  0.330 0.326 -1.04   0.0440  1.00   24973.   16489.
 5 beta[5]  -0.999  -1.00   0.828 0.822 -2.37   0.364   1.00   26126.   15510.
 6 beta[6]  -0.0637 -0.0680 0.845 0.854 -1.43   1.31    1.00   31265.   15003.
 7 beta[7]   1.54    1.55   0.664 0.663  0.449  2.61    1.00   27291.   15723.
 8 beta[8]   0.733   0.734  0.318 0.318  0.211  1.26    1.00   28048.   15306.
 9 beta[9]  -0.498  -0.485  0.792 0.802 -1.81   0.786   1.00   21930.   15612.
10 beta[10]  1.02    1.05   0.748 0.719 -0.265  2.20    1.00   15008.   13871.
# … with 38 more rows

Code

mfit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7989 0.8661 0.8909 0.8373 0.8476 0.8309 0.8887 0.8998

Trace plots

Code

mcmc_trace(mdrws["beta"])

Code

mcmc_trace(mdrws["gamma_site"])

Code

mcmc_trace(mdrws["gamma_epoch"])

Code

mcmc_trace(mdrws[c("tau_site", "tau_epoch")])

Posterior predictive

Code

y_ppc <- mdrws$y_ppc
ppc_dat <- bind_cols(acs_itt_nona_dat, tibble(y_ppc = y_ppc))

grp_ppc <- function(grp) {
  ppc_dat %>%
    group_by(grp = !!grp) %>%
    summarise(
      mean_y = mean(D28_death),
      rvar_mean_y_ppc = rvar_mean(y_ppc)
    )
}
plot_grp_ppc <- function(dat, lab = "") {
  dat %>%
    ggplot(aes(y = grp, xdist = rvar_mean_y_ppc)) +
    stat_interval(size = 2) +
    geom_point(aes(x = mean_y, y = 1:nrow(dat)), colour = "red", shape = 23) +
    labs(
      x = "Posterior predictive\nproportion", 
      y = lab)  
}

ppc_C <- grp_ppc(sym("CAssignment"))
ppc_ctry <- grp_ppc(sym("country"))
ppc_epoch <- grp_ppc(sym("epoch"))
ppc_site <- ppc_dat %>%
  group_by(Country = country, grp = site_raw) %>%
  summarise(
    mean_y = mean(D28_death),
    rvar_mean_y_ppc = rvar_mean(y_ppc)
  )

p1 <- plot_grp_ppc(ppc_C, "Anticoagulation\nintervention") + labs(x = "")
p2 <- plot_grp_ppc(ppc_ctry, "Country") + labs(x = "")
p3 <- plot_grp_ppc(ppc_epoch, "Epoch") + labs(x = "")
p4 <- plot_grp_ppc(ppc_site %>% filter(Country == "IN"), "Sites\nIndia")  + labs(x = "")
p5 <- plot_grp_ppc(ppc_site %>% filter(Country == "AU"), "Sites\nAustralia") + labs(x = "")
p6 <- plot_grp_ppc(ppc_site %>% filter(Country == "NP"), "Sites\nNepal")
p7 <- plot_grp_ppc(ppc_site %>% filter(Country == "NZ"), "Sites\nNZ")
p <- ((p3 | p1 / p2) / 
        ( (p4 | p5) / (p6 | p7) + 
            plot_layout(heights = c(3, 1)) ) ) +
  plot_layout(heights = c(1, 1.5), guides = "collect")
pth <- file.path("outputs", "figures", "outcomes", "secondary",
                 "7-3-primary-model-acs-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.5)
p

Posterior predictive distribution versus observed proportion (red diamond).

Sensitivity: Reduced model (no site terms)

The primary model is re-fit without the hierarchical site terms.

Fit model

mdat <- make_stan_data(
  acs_itt_nona_dat,
  c("inelgc3", "agegte60", "ctry"),
  "D28_death",
  c(10, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_epoch_only$sample(
    data = mdat,
    seed = 14646,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0,
    adapt_delta = 0.99
))
mdrws <- as_draws_rvars(mfit$draws(c("beta", "gamma_epoch", "tau_epoch")))
names(mdrws$beta) <- colnames(mdat$X)
epoch_map <- acs_itt_nona_dat %>% dplyr::count(epoch, epoch_lab)
names(mdrws$gamma_epoch) <- epoch_map$epoch_lab
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:4]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic")),
  setNames(exp(mdrws$beta[-(1:6)]), c("Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal"))
)

Odds ratio summary

odds_ratio_summary_table(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate	0.83	(0.42, 1.61)	0.88 (0.31)	0.71
Low with aspirin	0.84	(0.38, 1.83)	0.91 (0.37)	0.67
Therapeutic	3.39	(1.04, 10.68)	4.04 (2.61)	0.02
Ineligible aspirin	2.86	(0.85, 8.26)	3.30 (1.97)	0.05
Age ≥ 60	2.32	(1.28, 4.25)	2.44 (0.76)	0.00
Australia/New Zealand	0.76	(0.16, 3.05)	0.97 (0.78)	0.65
Nepal	4.56	(1.55, 13.75)	5.34 (3.29)	0.00

Odds ratio summary table for no epoch model.

Code

plot_epoch_terms(exp(mdrws$gamma_epoch))

Diagnostics

Code

mfit$summary(variables = c("beta", "gamma_epoch", "tau_epoch"))

# A tibble: 23 × 10
   variable    mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -5.01   -4.97   0.663 0.652 -6.16   -3.99    1.00   14715.   12970.
 2 beta[2]   1.01    1.00   0.478 0.478  0.224   1.78    1.00   20283.   15713.
 3 beta[3]  -0.246  -0.246  0.274 0.274 -0.701   0.201   1.00   22550.   15654.
 4 beta[4]  -0.548  -0.547  0.326 0.330 -1.08   -0.0154  1.00   19992.   15215.
 5 beta[5]  -0.936  -0.932  0.811 0.803 -2.30    0.386   1.00   25416.   16103.
 6 beta[6]  -0.0529 -0.0600 0.856 0.852 -1.45    1.37    1.00   26701.   13969.
 7 beta[7]   1.03    1.05   0.585 0.581  0.0344  1.96    1.00   25791.   13258.
 8 beta[8]   0.844   0.843  0.303 0.303  0.349   1.34    1.00   27646.   15597.
 9 beta[9]  -0.306  -0.278  0.756 0.753 -1.59    0.900   1.00   24776.   15176.
10 beta[10]  1.52    1.52   0.557 0.553  0.602   2.43    1.00   15983.   14684.
# … with 13 more rows

Code

mfit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8841 0.9230 0.8853 0.8829 0.9267 0.8939 0.8181 0.8834

Trace plots

Code

mcmc_trace(mdrws["beta"])

Code

mcmc_trace(mdrws["gamma_epoch"])

Code

mcmc_trace(mdrws["tau_epoch"])

Sensitivity: Reduced model (no epoch terms)

The primary model is re-fit without the epoch terms.

Fit model

mdat <- make_stan_data(
  acs_itt_nona_dat,
  c("inelgc3", "agegte60", "ctry"),
  "D28_death",
  c(10, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_site_only$sample(
    data = mdat,
    seed = 123513,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0,
    adapt_delta = 0.99
))

mdrws <- as_draws_rvars(mfit$draws(c("beta", "gamma_site", "tau_site")))
site_map <- acs_itt_nona_dat %>% dplyr::count(site_num, site)
site_facet <-   factor(mdat$region_by_site, labels = c("India", "Australia\nNew Zealand", "Nepal"))
names(mdrws$beta) <- colnames(mdat$X)
names(mdrws$gamma_site) <- site_map$site
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:4]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic")),
  setNames(exp(mdrws$beta[-(1:6)]), c("Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal"))
)

Odds ratio summary

odds_ratio_summary_table(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate	0.84	(0.42, 1.67)	0.89 (0.33)	0.69
Low with aspirin	0.88	(0.39, 1.95)	0.96 (0.41)	0.62
Therapeutic	2.93	(0.90, 9.11)	3.46 (2.22)	0.04
Ineligible aspirin	5.42	(1.37, 19.17)	6.63 (4.78)	0.01
Age ≥ 60	1.97	(1.07, 3.65)	2.07 (0.67)	0.01
Australia/New Zealand	0.59	(0.12, 2.64)	0.80 (0.71)	0.75
Nepal	2.28	(0.57, 7.75)	2.78 (2.03)	0.11

Odds ratio summary table for no epoch model.

Code

plot_site_terms(exp(mdrws$gamma_site), site_facet)

Diagnostics

Code

mfit$summary(variables = c("beta", "gamma_site", "tau_site"))

# A tibble: 39 × 10
   variable    mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -4.07   -4.05   0.519 0.504 -4.95   -3.26    1.00    8867.   12171.
 2 beta[2]   0.859   0.858  0.471 0.467  0.0896  1.63    1.00   22633.   14901.
 3 beta[3]  -0.221  -0.218  0.281 0.278 -0.685   0.240   1.00   26235.   16195.
 4 beta[4]  -0.505  -0.503  0.326 0.322 -1.05    0.0311  1.00   24775.   15954.
 5 beta[5]  -1.04   -1.04   0.834 0.841 -2.41    0.323   1.00   25838.   16057.
 6 beta[6]  -0.0811 -0.0795 0.849 0.844 -1.47    1.30    1.00   30912.   13477.
 7 beta[7]   1.67    1.69   0.667 0.654  0.568   2.75    1.00   23886.   15322.
 8 beta[8]   0.679   0.679  0.313 0.314  0.162   1.20    1.00   29941.   14381.
 9 beta[9]  -0.531  -0.521  0.790 0.786 -1.85    0.752   1.00   22361.   15767.
10 beta[10]  0.807   0.826  0.666 0.650 -0.323   1.86    1.00   13402.   14508.
# … with 29 more rows

Code

mfit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9495 0.8738 0.9165 0.8951 0.9422 0.8839 0.9512 0.8555

Trace plots

Code

mcmc_trace(mdrws["beta"])

Code

mcmc_trace(mdrws["gamma_site"])

Code

mcmc_trace(mdrws["tau_site"])

Sensitivity: Reduced model (no site or epoch)

The primary model is re-fit without the site and epoch terms.

Fit reduced model (no site or epoch)

mdat <- make_stan_data(
  acs_itt_nona_dat,
  c("inelgc3", "agegte60", "ctry"),
  "D28_death",
  c(10, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 12415,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0,
    adapt_delta = 0.99
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:4]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic")),
  setNames(exp(mdrws$beta[-(1:6)]), c("Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal"))
)

Odds ratio summary

odds_ratio_summary_table(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate	0.81	(0.41, 1.58)	0.86 (0.30)	0.73
Low with aspirin	1.26	(0.57, 2.68)	1.35 (0.55)	0.28
Therapeutic	2.58	(0.84, 7.33)	2.96 (1.73)	0.05
Ineligible aspirin	4.20	(1.28, 11.39)	4.75 (2.67)	0.01
Age ≥ 60	2.36	(1.30, 4.19)	2.46 (0.76)	0.00
Australia/New Zealand	0.57	(0.13, 2.18)	0.72 (0.55)	0.79
Nepal	1.97	(0.81, 4.39)	2.13 (0.93)	0.06

Odds ratio summary table for no epoch model.

Diagnostics

Code

mfit$summary(variables = c("beta"))

# A tibble: 10 × 10
   variable    mean  median    sd   mad      q5    q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -3.65   -3.65   0.257 0.252 -4.08   -3.23   1.00   11510.   13122.
 2 beta[2]   0.522   0.524  0.443 0.442 -0.215   1.24   1.00   13728.   15091.
 3 beta[3]  -0.274  -0.273  0.267 0.264 -0.711   0.167  1.00   16229.   13075.
 4 beta[4]  -0.631  -0.631  0.313 0.312 -1.14   -0.115  1.00   15687.   14173.
 5 beta[5]  -1.12   -1.11   0.794 0.783 -2.44    0.188  1.00   17669.   15548.
 6 beta[6]  -0.0715 -0.0679 0.831 0.815 -1.44    1.28   1.00   21121.   14415.
 7 beta[7]   1.41    1.43   0.559 0.552  0.451   2.28   1.00   18010.   14037.
 8 beta[8]   0.855   0.858  0.301 0.303  0.360   1.34   1.00   16662.   14081.
 9 beta[9]  -0.580  -0.566  0.725 0.729 -1.80    0.576  1.00   18319.   14119.
10 beta[10]  0.665   0.680  0.432 0.429 -0.0621  1.36   1.00   14109.   14057.

Code

mfit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9655 1.0383 0.9775 1.0953 1.0490 0.9577 1.0132 1.0193

Trace plots

Code

mcmc_trace(mdrws["beta"])

Sensitivity: Concurrent Controls

A sensitivity analysis is performed on the reduced concurrent control analysis sets for each intervention.

Intermediate

• Set: ACS-ITT-intermediate
• Covariates: anticoagulation intervention, age group, region of enrolment.

Overview of outcome

save_tex_table(
  make_mortality_table(acs_itt_concurc2_dat, "latex"),
  "outcomes/secondary/7-3-anticoagulation-concurrent-intermediate-summary")
make_mortality_table(acs_itt_concurc2_dat)

n (\%)	Patients	Known	Missing	Died by day 28
Low dose	610	596 (97.7)	14 (2.3)	19 (3.2)
Intermediate dose	613	603 (98.4)	10 (1.6)	15 (2.5)
Overall	1223	1199 (98.0)	24 (2.0)	34 (2.8)

Summary of outcome for intermediate dose concurrent enrolments.

Fit model

mdat <- make_stan_data_concurrent(
  acs_itt_concurc2_nona_dat,
  c("agegte60", "ctry"),
  "D28_death",
  c(2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 38135,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0,
    adapt_delta = 0.95
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]))
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate")),
  setNames(exp(mdrws$beta[-(1:2)]), c("Age \u2265 60", "Australia/New Zealand", "Nepal"))
)

Odds ratio summary

save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-3-primary-model-acs-itt-concurrent-intermediate-summary-table")
odds_ratio_summary_table(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate	0.78	(0.39, 1.54)	0.83 (0.29)	0.76
Age ≥ 60	2.84	(1.44, 5.62)	3.01 (1.08)	0.00
Australia/New Zealand	0.46	(0.12, 1.41)	0.54 (0.34)	0.91
Nepal	1.31	(0.46, 3.22)	1.45 (0.72)	0.30

Odds ratio summary table for model.

Posterior contrast

p <- plot_or_densities(mdrws$compare)
p

Diagnostics

Code

mfit$summary(variables = c("beta"))

# A tibble: 5 × 10
  variable   mean median    sd   mad     q5    q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -3.95  -3.94  0.258 0.256 -4.39  -3.54   1.00   11239.   11056.
2 beta[2]  -0.175 -0.173 0.246 0.246 -0.582  0.226  1.00   18118.   13824.
3 beta[3]   1.04   1.04  0.346 0.344  0.473  1.61   1.00   12024.   11313.
4 beta[4]  -0.804 -0.766 0.633 0.625 -1.89   0.169  1.00   17707.   12582.
5 beta[5]   0.250  0.269 0.497 0.490 -0.601  1.04   1.00   16895.   13100.

Code

mfit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9724 1.0293 1.0048 0.9900 1.0677 1.0275 1.0277 1.0140

Trace plots

Code

mcmc_trace(mdrws["beta"])

Low-dose with aspirin

• Set: ACS-ITT-aspirin
• Covariates: anticoagulation intervention, age group, region of enrolment.

Overview of outcome

save_tex_table(
  make_mortality_table(acs_itt_concurc3_dat, "latex"),
  "outcomes/secondary/7-3-anticoagulation-concurrent-stdaspirin-summary")
make_mortality_table(acs_itt_concurc3_dat)

n (\%)	Patients	Known	Missing	Died by day 28
Low dose	299	291 (97.3)	8 (2.7)	11 (3.8)
Intermediate dose	298	293 (98.3)	5 (1.7)	12 (4.1)
Low dose with aspirin	283	281 (99.3)	2 (0.7)	10 (3.6)
Overall	880	865 (98.3)	15 (1.7)	33 (3.8)

Summary of outcome for low-dose with aspirin concurrent enrolments.

Fit model

mdat <- make_stan_data_concurrent(
  acs_itt_concurc3_nona_dat,
  c("agegte60", "ctry"),
  "D28_death",
  c(2.5, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 13578,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:3]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Low with aspirin vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs low with aspirin" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin")),
  setNames(exp(mdrws$beta[-(1:3)]), c("Age \u2265 60", "Australia/New Zealand"))
)

Odds ratio summary

save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-3-primary-model-acs-itt-concurrent-stdaspirin-summary-table")
odds_ratio_summary_table(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate	1.09	(0.48, 2.48)	1.19 (0.53)	0.42
Low with aspirin	0.89	(0.38, 2.07)	0.97 (0.44)	0.61
Age ≥ 60	2.45	(1.22, 4.87)	2.60 (0.95)	0.01
Australia/New Zealand	0.41	(0.08, 1.58)	0.52 (0.40)	0.89

Odds ratio summary table for model.

Posterior contrast

p <- plot_or_densities(mdrws$compare)
p

Diagnostics

Code

mfit$summary(variables = c("beta"))

# A tibble: 5 × 10
  variable    mean  median    sd   mad     q5    q95  rhat ess_bulk ess_tail
  <chr>      <dbl>   <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -3.57   -3.56   0.246 0.244 -3.99  -3.18   1.00   16232.   14278.
2 beta[2]  -0.148  -0.145  0.301 0.301 -0.649  0.338  1.00   22800.   14161.
3 beta[3]   0.0127  0.0187 0.300 0.299 -0.493  0.495  1.00   22693.   13702.
4 beta[4]   0.894   0.895  0.355 0.355  0.311  1.47   1.00   16812.   14811.
5 beta[5]  -0.913  -0.881  0.752 0.746 -2.21   0.269  1.00   22255.   13232.

Code

mfit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.097 1.026 1.149 1.007 1.063 1.109 1.036 1.038

Trace plots

Code

mcmc_trace(mdrws["beta"])

Therapeutic

• Set: ACS-ITT-therapeutic
• Covariates: anticoagulation intervention, age group, region of enrolment.

Overview of outcome

save_tex_table(
  make_mortality_table(acs_itt_concurc4_dat, "latex"),
  "outcomes/secondary/7-3-anticoagulation-concurrent-therapeutic-summary")
make_mortality_table(acs_itt_concurc4_dat)

n (\%)	Patients	Known	Missing	Died by day 28
Low dose	79	75 (94.9)	4 (5.1)	3 (4.0)
Intermediate dose	65	63 (96.9)	2 (3.1)	1 (1.6)
Therapeutic dose	50	50 (100.0)	0 (0.0)	6 (12.0)
Overall	194	188 (96.9)	6 (3.1)	10 (5.3)

Summary of outcome for therapeutic dose concurrent enrolments.

Fit model

mdat <- make_stan_data_concurrent(
  acs_itt_concurc4_nona_dat,
  c("agegte60", "ctry"),
  "D28_death",
  c(2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 13578,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:3]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Therapeutic vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs therapeutic" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Therapeutic")),
  setNames(exp(mdrws$beta[-(1:3)]), c("Age \u2265 60", "India", "Australia/New Zealand"))
)

Odds ratio summary

save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-3-primary-model-acs-itt-concurrent-therapeutic-summary-table")
odds_ratio_summary_table(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate	0.66	(0.13, 3.11)	0.91 (0.86)	0.69
Therapeutic	2.64	(0.70, 10.68)	3.38 (2.75)	0.07
Age ≥ 60	2.80	(0.78, 11.41)	3.60 (3.05)	0.06
India	0.42	(0.08, 1.81)	0.55 (0.46)	0.87
Australia/New Zealand	0.28	(0.06, 1.02)	0.35 (0.26)	0.97

Odds ratio summary table for model.

Posterior contrast

p <- plot_or_densities(mdrws$compare)
p

Diagnostics

Code

mfit$summary(variables = c("beta"))

# A tibble: 6 × 10
  variable   mean median    sd   mad      q5    q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl> <dbl> <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -3.33  -3.29  0.566 0.551 -4.32   -2.47   1.00   12097.   10881.
2 beta[2]   0.986  0.969 0.555 0.543  0.0946  1.93   1.00   18013.   14124.
3 beta[3]  -0.228 -0.220 0.526 0.520 -1.10    0.622  1.00   19625.   14656.
4 beta[4]   1.04   1.03  0.676 0.660 -0.0390  2.17   1.00   12790.   10806.
5 beta[5]  -0.885 -0.861 0.784 0.779 -2.22    0.362  1.00   20168.   13677.
6 beta[6]  -1.29  -1.26  0.711 0.711 -2.50   -0.179  1.00   20226.   13969.

Code

mfit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0004 1.0536 1.0274 1.0999 1.0111 0.9799 1.0311 0.9898

Trace plots

Code

mcmc_trace(mdrws["beta"])